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Abstract 

In this paper we present efficient computational and symbolic algorithms for solving a nearly penta- 
diagonal linear systems. The implementation of the algorithms using Computer Algebra Systems (CAS) 
such as MAPLE, MACSYMA, MATHEMATICA, and MATLAB is straightforward. Two examples are 
given in order to illustrate the algorithms. 
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1 Introduction 

Many problems in mathematics and applied science require the solution of linear systems having nearly 
pentadiagonal coefficient matrices. This kind of linear system arises in many fields of numerical computa- 
tion and differential Equations [1, 2,3]. This article is a general case of the author article [3]. 

The main goal of the current paper is to develop an efficient algorithms for solving a general nearly 
pentadiagonal linear systems of the form: 

AX = Y 
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(1.2) 



X = (xi,x 2 , ■ ■ ■ ,x n ) T , Y = (yi,y 2 , ■ ■ ■ ,y n ) T and n > 5. 

A general nxn nearly pentadiagonal matrix A of the form (1.2) can be stored in bn — 4 memory locations by 
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using five vectors a = (0,1,02,... ,a n -2), a = (01, 02, . . . , a n -i,s), b = (t, b 2 , 63, ■ • • , b n ), b = (63, 64, ■ • ■ , b n ), 
and d = (di, d,2, ■ . ■ , d n ). When considering the system (1.1) it is advantageous to introduce three addi- 
tional vectors c = (a, C2, . . . , c„), e = (ei, e 2 , . . . , e„) and / = (fi,f 2 , fs, . . . , /„). These vectors are related 
to the vectors a, a, d, b, and 6. 

The current paper is organized as follows. In section 2, the main results are given. Illustrative examples 
are presented in section 3. In section 4, a conclusion is given. 



2 Main results 

In this section we are going to formulate a new computational and symbolic algorithms for solving a gen- 
eral nearly pentadiagonal linear systems of the form (1.1). To do this, we begin by considering the LU 
decomposition [4] of the matrix A in the form: 
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(2.1) 



From (2.1) we obtain 



r di 

d 2 - hei 

di — fiCi-i — -^-^di-2, a 2 = e r , 

. dn — flOn-2 — fn&n-l 



if i = 1 

if i = 2 

if i = 3(l)n-l 
if i = n, 



(2.2) 
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62 if i = 2 



j ' ' ' ' (bi-^)/ci-i if* = 3(l)n-l (2 ' 4) 

(&n ~ - /ie n -2)/c„-l if Z = 71 . 

It is not difficult to prove that the LU decomposition (2.1) exists only if Cj ^ 0, i = l(l)n — 1. Moreover, 
a general nearly pentadiagonal linear system (1.1) possesses a unique solution if, in addition, c„ ^ 0. On 
the other hand, the determinant of the matrix A is given by: 

n 

det(A) =J[ci, (2.5) 

i=l 

and this shows the importance of the vector c [5] . 
We may now formulate the following results. 

Algorithm 2.1 To solve the general nearly pentadiagonal linear system (1.1), we may proceed as follows: 
step 1: Set c\ — d\. 

step 2: If a = 0, then OUTPUT( ; the method is fails'); STOP, 
step 3: Set a n = s, 61 = t, e\ = a\, f 2 = j 2 -, and C2 = d 2 — f2&\- 
step 4: Ifc 2 = 0, then OUTPUT('the method is fails'); STOP, 
step 5: Set e 2 = a 2 - /2ai, e„ = d 2 - fia n , and 0,2 = en- 
step 6: For i = 3, 4, . . . , n — 1 Compute 
fi = ipi - ^e i _ 2 )/c i _i ) 

^ = c?i - /»ei_i - ^-^-2, 

//Ci = 0, tfien OUTPUT('the method is fails'); STOP, 
Ifi = 3, then e { = a { - /jOj-i - efee e; = Oj - /;aj_i. 

step 7: Compute 

fi = (bn - ■sf^e„„ 3 )/c„_ 2 , 

,/n = (&n - — 3^n-3 - fl^n-l)/c n -l, 
Cn <^n fl&n—l fnCn — 1, 

step 8: Set z\ = y 1 , z 2 = yi - fiZ\. 
step 9: For i = 3, 4, . . . , n — 1 Compute 

Zi =Vi~ JiZi-i - 

step 10: Set z n = y n - f n z n -i - f\z n -i - -^—^n-z- 
step 11: Compute the solution vector x using 

t — Sn. v . — z„-l-e„-ia:„ 

For i = n — 2,n — 3, ... ,2 compute 

_ Zi-eiXi +1 -aiX i+ 2 

5et x\ = z i- e i x 2-ai x 3- sx 4 



The new algorithm 2.1 will be referred to as KNPENTA algorithm. KNPENTA algorithm for solving 
the nearly pentadiagonal system (1.1) is generally preferable because the conditions Cj 7^ 0, i = l(l)n are 
sufficient for its validity. The advantage of the vector c is now clear. 

The following symbolic algorithm is developed in order to remove the cases where the numeric algorithm 
KNPENTA fails. 

Algorithm 2.2 To solve the general nearly pentadiagonal linear system (1.1), we may proceed as follows: 
step 1: Set c\ — d\. 

step 2: Set c\ = x and d\ = x(x is just a symbolic name) whenever c\ = 0. 
step 3: Set a n = s, b\ = t, e\ = a\, f 2 = t 2 -, and c 2 = d 2 — f 2 e\. 
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step 4 : Set C2 = x whenever C2 = 0. 

step 5: Set e 2 = a2 — feai, e n = 5,2 — jia n , and 0,2 = e r , 

step 6: For i = 3, 4, . . . , n — 1 Compute 



bi ~ 



Ci = di - fai-i - -^ai-2, 
Set Ci — x whenever Ci — 0, 

Ifi = 3, then a = a { - f t a i - 1 - ^a„ else e t = a { - Uo-i-\- 
step 7: Compute 

fi = (bn - ■^e„_ 3 )/c„_ 2 , 

fn = (bn - ~ flC n ~2) / C n -\ , 

Cn d n f\(l n — 2 fnC-n—l; 

Set c n = x whenever c n = 0. 
step 8: Set z\ = yi, z 2 = y 2 - h%\- 
step 9: For i = 3, 4, . . . , n — 1 Compute 

Z%=Vi — fiZi-l - -^7 2 z i-2- 
step 10: Set z n = y n - f n z n -i - hz n -2 - 
step 11: Compute the solution vector x using 

For i = n — 2,n — 3, ... ,2 compute 

_ Zi-eiXi + i~aiX i+ 2 
x i r . 7 

Set X-\ = z i~ e i ;r 2-5ia3-sa:4 

step 12: Substitute x — in all expressions of the solution vector Xi,i — 1,2, ... ,n. 



The symbolic algorithm 2.2 will be referred to as KSNPENTA algorithm. 
In [6], Claerbout showed that the two-dimensional Laplacian operator, which appears in 3-D finite-difference 
migration, has the form of pentadiagonal matrix. If we choose di = —4, i = l(l)n, s = t = and 
at = bi = hi = bi = 1 Vi, we can obtain it. 



3 Illustrative Examples 

In this section we are going to give illustrative examples 

Example 3.1. Solve the nearly pentadiagonal linear system of size 10 given by 
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by using the KNPENTA algorithm and KSNPENTA algorithm. 



Solution 



(3-1) 



(i) The application of the KNPENTA algorithm gives: 



• ci = 3(Step 1). 
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• a n = 5, bi = -2, ei = — 1, / 2 = -J-, and c 2 = |(Step 3). 



• e 2 = 3, e„ = ^ and a 2 = -y (Step 5). 



-253 -105 4012 -368 2045671 



[/3> /4> /$» /e> hif%i /s>] I 4 ' ^' a*> J ' 22TI ' Mil • l,sc>l:«J' 

[„ „ „ „ „ „ ^ 1 — [35 1 552 757 -1313 14341 4T122S2 i 

[l-3, c 4, c 5, c 6, w, c 8, C 9J L 4 , 35 , 92 ' 1514 ' 303 ' 186433-1 

r„ _ _ _ _ _ c 1 — no O 934 1393 26873 -1371 -35.12041 ug+„„ c\ 

[e 3 ,e4,e5,e 6 ,e7,e 8 ,e 9 J - [12, 2, -gg-, - T§4 -, -3^-, -^j-, 186433 ](btcp b). 



f _ -91736 f _ -1203361 j _ 701215 (o t ^ 
•' 1 — 186433 i rn- — 4 n 2 262 ™ ^ " 19866 ' )■ 



• [21,22] = [30,33] (Step 8). 

.T^^^^^^^l r317 on 8609 11475 238883 138311 129996 1 fe+„„ Q\ 

• [Z3,24, 25, 26, 27, 28, 29j - [— ,32, -35-, - [Er , ^543-, , ^434T-J l bte P y )- 

• «i0 = niir ( Stc P 10 )- 

• [xi,x 2 ,x 3 ,X4,x 5 , x 6 , x 7 , x 8 , x 9 , x w ] = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] (Step 11). 

Also the determinant of the matrix A is det(A) = —145151505 by using (2.5). 

(ii) The application of the KSNPENTA algorithm gives: 

X := nearly _penta(b, b, d, a, a, y) — [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. 

Example 3.2. Solve the nearly pentadiagonal linear system of size 10 given by 
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by using the KNPENTA algorithm and KSNPENTA algorithm. 
Solution 

(i) The application of the KNPENTA algorithm gives: 

The method is broken down since c\ = d\ = 0. 

(ii) The application of the KSBPENTA algorithm gives: 

X := nearly_penta{b,b,d,a,a,y)— i L "' J2 ' JN ' 



(3.2) 



1490963 s-2728658 11371544 2,-12278961 



U589918x-4092987 ■ 1589918 2,-4092987 ' 4589918 2,-4092987 ' 
g 3279301 2,-2728658 o 7767312 2,-6821645 t /o 90274465 x-73673766 



- 3 • 



- 1 /3 ■ 
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i in 95213875 2,-85952727 i /r 1888517712,-163719480 
i° 4589918 2,-4092987 > V d ' 



1/15 



4589918 j— 1.092987 ' 
612846296 2,-552553245 



4589918 2,-4092987 



2/15 



342051290 .r - 300974025 1 



J2, = 



4589918 2,-4092987 ' "I 4589918 2,-4092987 

= [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. 
Also the determinant of the matrix A is det(A) = 61394805 and for more details about how to call this 
procedure, see appendix 1. 
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4 Conclusion 



The methods described here are very effective, provided that optimal LU factorization is used. Our sym- 
bolic algorithm is competitive with the other methods for solving a nearly pentadiagonal linear system 
which appears in many applications. 
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Appendix 1. A Maple procedure for solving a nearly pentadiagonal linear system 

> # A Maple Procedure. 

> # Written by Dr. A. A. Karawia 18-6-2008. 

> # To compute the solution of A general nearly pentadiagonal Linear system. 

> restart: 

> with(linalg,vector,vectdim): 

> nearly _penta:=proc(bb::vector,b::vector, d::vector,a::vector,aa::vector,y::vector) 
local i,j,k,n; global T,e,c,f,z,X; 

n:=vectdim(d):e:=array(l..n):c:=array(l..n): f:=array(l..n):z:=array(l..n):X:= array (l..n): 
^components of the vectors e, c, and f # 

c[l]:=d[l]:if c[l]=0 then c[l]:=x; d[l]:=x;fi: c[l]:=a[l]:f[2]:=simplify(b[2]/c[l]): 
c[2]:=simplify(d[2]-e[l]*f[2]):if c[2]=0 then c[2]:=x;fi: 
e[2]:=simplify(a[2]-f[2]*aa[l]):e[n]:=simplify(aa[2]-f[2]*a[n]): aa[2]:=e[n]: 
for i from 3 to n-1 do 

f [i] : ^simplify ( (b [i]-bb [i] *e [i-2] /c [i-2] ) /c [i- 1] ) : 

if i=3 then e[i]:=simplify(a[i]-f[i]*aa[i-l]-bb[i]*a[n]/c[i-2]); else e[i]:=simplify(a[i]-f[i]*aa[i-l]); 
fi: 

c[i]:=simplify(d[i]-bb[i]*aa[i-2]/c[i-2]-e[i-l]*f[i]); if c[i]=0 then c[i]:=x; fi: 
end do: 

f[l]:=simpIify((bb[n]-b[l]*e[n-3]/c[n-3])/c[n-2]): 
f[n]:=simplify((b[n]-b[l]*aa[n-3]/c[n-3]-f[l]*e[n-2])/c[n-l]): 
c[n]:=simplify(d[n]-f[l]*aa[n-2]-c[n-l]*f[n]):if c[n]=0 then c[n]:=x; fi: 

# To compute the vector Z # 
z[l]:=y[l]:z[2]:=y[2]-f[2]*z[l]:i:='i': 

for i from 3 to n-1 do 

z[i]:=simplify(y[i]-bb[i]*z[i-2]/c[i-2]-f[i]*z[i-l]): 

end do: 

z[n]:=simplify(y[n]-f[n]*z[n-l]-f[l]*z[n-2]-b[l]*z[n-3]/c[n-3]): 

# To compute the Solution of the system X. # 
X[n]:=z[n]/c[n]:i:='i': 

X[n-l]:=simplify((z[n-l]-e[n-l]*X[n])/c[n-l]): 
for i from n-2 by -1 to 2 do 

X[i]:=simplify((z[i]-e[i]*X[i+l]-aa[i]*X[i+2])/c[i]): 

end do: 

X[l]:=simplify((z[l]-e[l]*X[2]-aa[l]*X[3]-a[n]*X[4])/c[i]): 

=ff To compute the determinant T ^ 

T:=subs(x=0, simplify (product (c[r],r=l..n))); 

eval(X): 

end: 

> # Call no. 1 for the procedure nearly _pcnta. # 

> x:='x': 

> aa:=aa:=vector([3,2, 1,3, 1,-5,-4, 20]); 
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aa := [3, 2, 1, 3, 1, -5, -4, 20] 
>a:=vector([-l, 1,5,1,5,7,3,1,3,5]); 

a:=[-l, 1,5,1,5,7,3,1,3,5] 

> d:=vector([3, 2,5, 1,2, 2, 12,3,21,31]); 

d:=[3,2, 5, 1,2,2, 12,3, 21,31] 

> b:=vector([-2,-2,-4,-2,l,-3,l,5,ll,-9]); 

b:=[-2,-2,-4,-2,l,-3,l,5,ll,-9] 

> bb:=vector([0,0,3,3,6,3,-8, 2,3,4]); 

bb:=[0,0,3,3,6,3,-8,2,3,4] 

> y:=vector([30, 13,35, 27,69, 18,38,280,328,247]); 

y:=[30,13,35,27,69, 18,38,280,328,247] 

> X:=nearly_penta(bb,b,d,a,aa,y); 

X := [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] 

>T; 

-145151505 

> # End of call no. 1. # 



> # Call no. 2 for the procedure nearly _penta. # 

> x:='x': 

> aa:=vector([3, 2, 1,3, 1,-5,-4,20]); 

aa := [3, 2, 1, 3, 1, -5, -4, 20] 
>a:=vector([-l,l,5,l,5,7,3,l,3,5]); 

a:=[-l,l,5,l,5,7,3,l,3,5] 

> d:=vector([0, 2,5, 1,2, 2, 12,3,21,31]); 

d:=[0,2, 5, 1,2,2, 12,3, 21,31] 

> b:=vector([-2,-2,-4,-2, 1,-3, 1,5, 11,-9]); 

b:=[-2,-2,-4,-2,l,-3,l,5,ll,-9] 

> bb:=vector([0,0,3,3,6,3,-8, 2,3,4]); 

bb:=[0,0,3,3,6,3,-8,2,3,4] 

> y:=vector([27,13,35, 27,69, 18,38,280,328,247]); 

y:=[27,13,35,27,69, 18,38,280,328,247] 

> X:=nearly_penta(bb,b,d,a,aa,y); 

Y ._ [ -4092987 o 1490963 rr-2728658 11371544 z-12278961 c 3279301 z-2728658 o 7767312 z-6821645 

U589918X-4092987 ' ° 4589918 x-4092987' 4589918 z-4092987 ' 4589918 z-4092987 ' ° 4589918 z-4092987 ' 
i /q 90274465 a-73673766 i /q 95213875^-85952727 i /c 188851771 z-163719480 i /i c 612846296 z-552553245 

'° 4589918 x-4092987 ' ' ° 4589918 2-4092987 ' ' ° 4589918 s-4092987 ' / ° 4589918 2-4092987 ' 
o/i c 342051296 rr-306974025 1 
A l 10 4589918 rr-4092987 -I 

> T; 

61394805 

> x:=0:X:=map(eval,op(X)); 

X := [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] 

> # End of call no. 2. # 
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